Exhausted CD8 T cells are associated with ICB responsiveness in breast cancer
Background
The single-cell dataset generated by Bassez et al. (2021) Nature Medicine contains single-cell data from breast carcinoma biopsies from 29 patients taken immediately before (pre-treatment) and 9 +/- 2 days after (on-treatment) receiving anti-PD1 immunotherapy. Lacking an objective measure of clinical benefit following therapy, in this study the authors use T cell expansion following treatment as a proxy for response.
Goal
To compare the T cell subtype composition between responders and non-responders using ProjecTILs and a reference TIL atlas
R Environment
# renv::activate()
renv::restore()
# remotes::install_github('carmonalab/ProjecTILs', ref='1.0.0')
# remotes::install_version('Seurat', version = '4.0.1')
library(ProjecTILs)
library(patchwork)
library(ggplot2)
library(reshape2)
library(parallelly)
options(parallelly.makeNodePSOCK.setup_strategy = "sequential") # workaround to issue https://github.com/HenrikBengtsson/future/issues/511scRNA-seq data preparation
ddir <- "input/Bassez_breast_data"
if (!file.exists(ddir)) {
dir.create(ddir)
dataUrl <- "https://drive.switch.ch/index.php/s/7jdqnvPZuJrriZB/download"
download.file(dataUrl, paste0(ddir, "/tmp.zip"))
unzip(paste0(ddir, "/tmp.zip"), exdir = ddir)
file.remove(paste0(ddir, "/tmp.zip"))
}
# Count matrices
f1 <- sprintf("%s/1864-counts_tcell_cohort1.rds", ddir)
cohort1 <- readRDS(f1)
dim(cohort1)[1] 25288 53382# Metadata
meta1 <- read.csv(sprintf("%s/1870-BIOKEY_metaData_tcells_cohort1_web.csv", ddir))
rownames(meta1) <- meta1$Cell
head(meta1) Cell nCount_RNA
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1 BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1 1252
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1 BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1 1869
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1 BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1 1000
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1 BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1 1288
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1 BIOKEY_13_Pre_AAAGATGGTTTGACAC-1 2056
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1 BIOKEY_13_Pre_AAAGATGTCAACGCTA-1 1224
nFeature_RNA patient_id timepoint expansion
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1 700 BIOKEY_13 Pre n/a
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1 995 BIOKEY_13 Pre n/a
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1 627 BIOKEY_13 Pre n/a
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1 681 BIOKEY_13 Pre n/a
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1 789 BIOKEY_13 Pre n/a
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1 634 BIOKEY_13 Pre n/a
BC_type cellSubType cohort
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1 HER2+ gdT treatment_naive
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1 HER2+ NK_REST treatment_naive
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1 HER2+ CD4_N treatment_naive
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1 HER2+ CD8_EM treatment_naive
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1 HER2+ CD8_N treatment_naive
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1 HER2+ CD8_RM treatment_naivedata.seurat <- CreateSeuratObject(cohort1, project = "Cohort1_IT", meta.data = meta1)Subset pre-treatment biopsies
table(data.seurat$timepoint)
On Pre
27528 25854 data.seurat <- subset(data.seurat, subset = timepoint == "Pre")Subset T cells according to annotation by the authors
table(data.seurat$cellSubType)
CD4_EM CD4_EX CD4_EX_Proliferating
4810 1707 107
CD4_N CD4_REG CD4_REG_Proliferating
3007 2861 69
CD8_EM CD8_EMRA CD8_EX
5443 290 2286
CD8_EX_Proliferating CD8_N CD8_RM
340 354 1963
gdT NK_CYTO NK_REST
1000 230 938
Vg9Vd2_gdT
449 data.seurat <- subset(data.seurat, subset = cellSubType %in% c("NK_REST", "Vg9Vd2_gdT",
"gdT", "NK_CYTO"), invert = T)Basic QC, remove outliers in # counts, # detected genes, and mitochondrial and ribosomal content
percent.ribo.dv <- PercentageFeatureSet(data.seurat, pattern = "^RP[LS]")
percent.mito.dv <- PercentageFeatureSet(data.seurat, pattern = "^MT-")
data.seurat <- AddMetaData(data.seurat, metadata = percent.ribo.dv, col.name = "percent.ribo")
data.seurat <- AddMetaData(data.seurat, metadata = percent.mito.dv, col.name = "percent.mito")
Idents(data.seurat) <- "patient_id"
VlnPlot(data.seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.ribo", "percent.mito"),
ncol = 2, pt.size = 0)
data.seurat <- subset(data.seurat, subset = nFeature_RNA > 500 & nFeature_RNA < 6000 &
nCount_RNA > 600 & nCount_RNA < 25000 & percent.ribo < 50 & percent.mito < 15)Remove samples that are too small and downsample large ones, in order to have similar contribution from all patients/samples contribute. Downsampling large samples also speeds up downstream computations.
minCells <- 500
ds <- 1000
Idents(data.seurat) <- "patient_id"
table(data.seurat$patient_id)
BIOKEY_1 BIOKEY_10 BIOKEY_11 BIOKEY_12 BIOKEY_13 BIOKEY_14 BIOKEY_15 BIOKEY_16
3014 1006 543 2511 1220 920 1026 2070
BIOKEY_17 BIOKEY_18 BIOKEY_19 BIOKEY_2 BIOKEY_20 BIOKEY_21 BIOKEY_22 BIOKEY_23
96 182 1430 945 109 194 68 35
BIOKEY_24 BIOKEY_25 BIOKEY_26 BIOKEY_27 BIOKEY_28 BIOKEY_29 BIOKEY_3 BIOKEY_30
339 95 85 594 575 162 264 61
BIOKEY_31 BIOKEY_4 BIOKEY_5 BIOKEY_6 BIOKEY_7 BIOKEY_8 BIOKEY_9
346 2282 1688 606 63 135 168 largeSamples <- names(table(data.seurat$patient_id)[table(data.seurat$patient_id) >=
minCells])
data.seurat <- subset(data.seurat, idents = largeSamples)
table(data.seurat$patient_id)
BIOKEY_1 BIOKEY_10 BIOKEY_11 BIOKEY_12 BIOKEY_13 BIOKEY_14 BIOKEY_15 BIOKEY_16
3014 1006 543 2511 1220 920 1026 2070
BIOKEY_19 BIOKEY_2 BIOKEY_27 BIOKEY_28 BIOKEY_4 BIOKEY_5 BIOKEY_6
1430 945 594 575 2282 1688 606 data.seurat <- subset(data.seurat, cells = WhichCells(data.seurat, downsample = ds))
table(data.seurat$patient_id)
BIOKEY_1 BIOKEY_10 BIOKEY_11 BIOKEY_12 BIOKEY_13 BIOKEY_14 BIOKEY_15 BIOKEY_16
1000 1000 543 1000 1000 920 1000 1000
BIOKEY_19 BIOKEY_2 BIOKEY_27 BIOKEY_28 BIOKEY_4 BIOKEY_5 BIOKEY_6
1000 945 594 575 1000 1000 606 ProjecTILs analysis
Here we will project the query data onto a reference TIL atlas to interpret cell states in the context of reference T cell subtypes.
For optimal batch-effect correction, we prefer projecting each patient/batch separately (although results tend to be very similar if projecting all together)
Here we use filter.cells=F to save processing time, because we rely on the authors’ cell type annotation to subset T cells.
For large samples make.projection can take several minutes. Set the number of parallel jobs with ncores according to your computational resources (you can check future::availableCores(), although parallelization is usually limited by available RAM)
ncores = 5
ref <- load.reference.map() # by default, murine reference TIL atlas v1.0[1] "Loading Default Reference Atlas..."
[1] "/Users/admin/gitHub/ProjecTILs_CaseStudies/ref_TILAtlas_mouse_v1.rds"
[1] "Loaded Reference map ref_TILAtlas_mouse_v1"data.split <- SplitObject(data.seurat, split.by = "patient_id")
query.projected <- make.projection(data.split, ref, filter.cells = F, ncores = ncores)[1] "Using assay RNA for BIOKEY_13"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_13 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_10"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_10 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_16"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_16 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_14"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_14 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_19"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_19 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_28"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_28 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_15"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_15 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_5"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_5 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_12"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_12 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_1"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_1 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_4"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_4 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_11"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_11 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_2"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_2 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_6"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_6 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_27"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_27 to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"See query cells for two selected the samples, projected onto the reference atlas
plot.projection(ref, query = query.projected$BIOKEY_13) + ggtitle("Projection sample 1")
plot.projection(ref, query = query.projected$BIOKEY_10) + ggtitle("Projection sample 2")
For simplified downstream analysis, merge back ProjecTILs results list of objects into a single Seurat object
query.projected <- lapply(query.projected, function(x) {
cellstate.predict(ref = ref, query = x, reduction = "umap", ndim = 2)
})
query.projected.merged <- suppressMessages(Reduce(ProjecTILs:::merge.Seurat.embeddings,
query.projected))We visually inspect now how the gene profiles of projected cells match those of the reference, for each T cell subtype.
A query dataset might contain a subtype that is not well represented in the atlas. That would requiere splitting heterogeneous populations after projection by higher-dimensions (e.g. Find Discriminant Dimensions ) In this case, radar/spider profiles look fine, meaning that in general terms, current reference atlas is capturing the main features of the query data:
plot.states.radar(ref, query = query.projected.merged, min.cells = 30, genes4radar = c("Foxp3",
"Cd4", "Cd8a", "Tcf7", "Ccr7", "Sell", "Il7r", "Gzmb", "Gzmk", "Pdcd1", "Lag3",
"Havcr2", "Tox", "Entpd1", "Cxcr3", "Cxcr5", "Ifng", "Cd200", "Xcl1", "Itgae")) #,'Cxcl13',
We can also verify that gene profiles are consistent across samples/patients, i.e. batch effects have been mitigated
plot.states.radar(ref, query = query.projected[1:5], min.cells = 30, genes4radar = c("Foxp3",
"Cd4", "Cd8a", "Tcf7", "Ccr7", "Sell", "Il7r", "Gzmb", "Gzmk", "Pdcd1", "Lag3",
"Havcr2", "Tox", "Entpd1", "Cxcr3", "Ifng", "Cd200", "Xcl1", "Itgae"))
Compare responders vs. non-responders (based on clonal expansion) - pre-therapy
# Split by expansion (defined by authors)
query.sub.byExp <- subset(query.projected.merged, subset = expansion %in% c("E",
"NE"))
query.sub.byExp <- SplitObject(query.sub.byExp, split.by = "expansion")
n_e <- length(table(query.sub.byExp$E$patient_id))
n_ne <- length(table(query.sub.byExp$NE$patient_id))
a <- plot.projection(ref, query = query.sub.byExp$E, linesize = 0.5, pointsize = 0.2) +
ggtitle(paste0("Responders; n=", n_e))
b <- plot.projection(ref, query = query.sub.byExp$NE, linesize = 0.5, pointsize = 0.2) +
ggtitle(paste0("Non-Responders; n=", n_ne))
a
b
c <- plot.statepred.composition(ref, query.sub.byExp$E, metric = "percent") + ylim(0,
30) + ggtitle(paste0("Responders; n=", n_e))
d <- plot.statepred.composition(ref, query.sub.byExp$NE, metric = "percent") + ylim(0,
30) + ggtitle(paste0("Non-Responders; n=", n_ne))
c | d
The analysis shows that ICB-responding tumors have a higher proportion of exhausted (Tex) and precursor-exhausted (Tpex) CD8 TILs, as well as a lower ratio of Treg/Tex.
Calculate and visualize fold-changes of TIL composition in responders vs non-responders
which.types <- table(query.projected.merged$functional.cluster) > 20 # disregard very small populations
stateColors_func <- c("#edbe2a", "#A58AFF", "#53B400", "#F8766D", "#00B6EB", "#d1cfcc",
"#FF0000", "#87f6a5", "#e812dd")
states_all <- levels(ref$functional.cluster)
names(stateColors_func) <- states_all
cols_use <- stateColors_func[names(which.types)][which.types]
query.projected.merged$functional.cluster <- factor(query.projected.merged$functional.cluster,
levels = states_all)
query.list <- SplitObject(query.projected.merged, split.by = "expansion")
norm.c <- table(query.list[["NE"]]$functional.cluster)/sum(table(query.list[["NE"]]$functional.cluster))
norm.q <- table(query.list[["E"]]$functional.cluster)/sum(table(query.list[["E"]]$functional.cluster))
foldchange <- norm.q[which.types]/norm.c[which.types]
foldchange <- sort(foldchange, decreasing = T)
tb.m <- melt(foldchange)
colnames(tb.m) <- c("Cell_state", "Fold_change")
pll <- list()
ggplot(tb.m, aes(x = Cell_state, y = Fold_change, fill = Cell_state)) + geom_bar(stat = "identity") +
scale_fill_manual(values = cols_use) + geom_hline(yintercept = 1) + scale_y_continuous(trans = "log2") +
theme(axis.text.x = element_blank(), legend.position = "left") + ggtitle("Responder vs. Non-responder")
Tumors in ICB-responders have a ~4 fold-change increase in Tex and a ~2 fold-change increase in both Tpex and Tregs.
Discussion
Our automated ProjecTILs analysis of this patient cohort indicate that compared to non-responders, early responders to PD-1 blockade (in terms of T cell expansion) have higher relative amounts of CD8 Tex and Tregs infiltrating their tumors, and a lower Treg to T(p)ex ratio, in agreement with previous studies (Daud et al. 2016; Thommen et al. 2018; Kumagai et al. 2020; Bassez et al. 2020). This supports the notion that patients with a pre-existing tumor-specific CD8 T cell response in their tumors are more likely to respond to PD-1 blockade therapy.
Because tumoricidal CD8 Tex are derived from Tpex (Siddiqui et al. 2019; Miller et al. 2019), the presence of CD8 Tex seems to be a good proxy for the presence of smaller populations of quiescent Tpex, that have the ability to proliferate and differentiate into Tex cells following PD-1 blockade (Siddiqui et al. 2019; Miller et al. 2019). Indeed we observed a small population of Tpex-like cells, with increased levels R vs RN. Because there are very few Tpex their gene profile might be less robust than others, but compared to Tex they display higher levels of Tpex-associated genes, such as XCL1 and CD200, and lower levels of terminal-exhaustion markers such as HAVCR2 and ENTPD1.
In their study, Bassez et al. also observed increased levels in R vs NR of a population they referred to as exhausted CD4 T cells. Here we observed also a trend for increased levels of follicular-helper CD4 T cells, that display exhaustion features (e.g. TOX, PDCD1, ENTPD1 expression), although the most prominent enrichment in ICB-responders is for exhausted CD8 T cells.
Further reading
Dataset original publication - Bassez et al
ProjecTILs case studies - INDEX - Repository
The ProjecTILs method Andreatta et. al (2021) Nat. Comm. and code
References
Bassez, A., Vos, H., Van Dyck, L. et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med 27, 820–832 (2021). https://doi.org/10.1038/s41591-021-01323-8
Daud, A.I., Loo, K., Pauli, M.L., Sanchez-Rodriguez, R., Sandoval, P.M., Taravati, K., Tsai, K., Nosrati, A., Nardo, L., Alvarado, M.D., Algazi, A.P., Pampaloni, M.H., Lobach, I.V., Hwang, J., Pierce, R.H., Gratz, I.K., Krummel, M.F., Rosenblum, M.D., 2016. Tumor immune profiling predicts response to anti-PD-1 therapy in human melanoma. J. Clin. Invest. 126, 3447–3452. https://doi.org/10.1172/JCI87324
Gros, A., Robbins, P.F., Yao, X., Li, Y.F., Turcotte, S., Tran, E., Wunderlich, J.R., Mixon, A., Farid, S., Dudley, M.E., Hanada, K., Almeida, J.R., Darko, S., Douek, D.C., Yang, J.C., Rosenberg, S. a, 2014. PD-1 identifies the patient-specific in filtrating human tumors. J. Clin. Invest. 124, 2246–59. https://doi.org/10.1172/JCI73639.2246
Miller, B.C., Sen, D.R., Al Abosy, R., Bi, K., Virkud, Y.V., LaFleur, M.W., Yates, K.B., Lako, A., Felt, K., Naik, G.S., Manos, M., Gjini, E., Kuchroo, J.R., Ishizuka, J.J., Collier, J.L., Griffin, G.K., Maleri, S., Comstock, D.E., Weiss, S.A., Brown, F.D., Panda, A., Zimmer, M.D., Manguso, R.T., Hodi, F.S., Rodig, S.J., Sharpe, A.H., Haining, W.N., 2019. Subsets of exhausted CD8+ T cells differentially mediate tumor control and respond to checkpoint blockade. Nat. Immunol. 20, 326–336. https://doi.org/10.1038/s41590-019-0312-6
Siddiqui, I., Schaeuble, K., Chennupati, V., Fuertes Marraco, S.A., Calderon-Copete, S., Pais Ferreira, D., Carmona, S.J., Scarpellino, L., Gfeller, D., Pradervand, S., Luther, S.A., Speiser, D.E., Held, W., 2019. Intratumoral Tcf1+PD-1+CD8+ T Cells with Stem-like Properties Promote Tumor Control in Response to Vaccination and Checkpoint Blockade Immunotherapy. Immunity 50, 195–211.e10. https://doi.org/10.1016/J.IMMUNI.2018.12.021
Thommen, D.S., Koelzer, V.H., Herzig, P., Roller, A., Trefny, M., Dimeloe, S., Kiialainen, A., Hanhart, J., Schill, C., Hess, C., Prince, S.S., Wiese, M., Lardinois, D., Ho, P.C., Klein, C., Karanikas, V., Mertz, K.D., Schumacher, T.N., Zippelius, A., 2018. A transcriptionally and functionally distinct pd-1 + cd8 + t cell pool with predictive potential in non-small-cell lung cancer treated with pd-1 blockade. Nat. Med. 24, 994. https://doi.org/10.1038/s41591-018-0057-z
Kumagai, S., Togashi, Y., Kamada, T. et al. The PD-1 expression balance between effector and regulatory T cells predicts the clinical efficacy of PD-1 blockade therapies. Nat Immunol 21, 1346–1358 (2020). https://doi.org/10.1038/s41590-020-0769-3